Surface emitter

In [1]:

# Useful for debugging
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'

import os
import yaml
import numpy as np

# Useful for debugging %load_ext autoreload %autoreload 2 %config InlineBackend.figure_format = 'retina' import os import yaml import numpy as np

In [2]:

from distgen import Generator

from distgen import Generator

In [3]:

input_yaml = """
n_particle: 100000
species: electron
start: 
  type: free
random:
  type: hammersley
total_charge:
  units: pC
  value: 10
r_dist:
  max_r:
    units: mm
    value: 1
  type: radial_uniform
p_dist:
  sigma_p:
    value: 1
    units: eV/c
  avg_p: 
    value: 100
    units: eV/c
  type: gaussian
"""

input_yaml = """ n_particle: 100000 species: electron start: type: free random: type: hammersley total_charge: units: pC value: 10 r_dist: max_r: units: mm value: 1 type: radial_uniform p_dist: sigma_p: value: 1 units: eV/c avg_p: value: 100 units: eV/c type: gaussian """

In [4]:

gen = Generator(input_yaml, verbose=True)

gen = Generator(input_yaml, verbose=True)

In [5]:

P = gen.run()

P = gen.run()

Distribution format: None
   Warning: no output file specified, defaulting to "None".
Output file: None

Creating beam distribution....
   Beam starting from: free
   Total charge: 10 pC.
   Number of macroparticles: 100000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 1 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   p distribution: Gaussian
      avg_p = 100 eV/c, sigma_p = 1.000 eV/c
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   Shifting avg_x = -4.85224E-08 mm -> 0 mm
   Scaling sigma_x = 0.499991 mm -> 0.5 mm
   Shifting avg_y = 1.44675E-06 mm -> 0 mm
   Scaling sigma_y = 0.49999 mm -> 0.5 mm
   Shifting avg_px = -3.3074E-05 eV/c -> -6.12323E-15 eV/c
   Scaling sigma_px = 57.7382 eV/c -> 57.7379 eV/c
   Shifting avg_py = 1.54491E-06 eV/c -> 0 eV/c
   Scaling sigma_py = 57.7382 eV/c -> 57.7379 eV/c
   Shifting avg_pz = -0.000118434 eV/c -> 6.12323E-15 eV/c
   Scaling sigma_pz = 57.7372 eV/c -> 57.7379 eV/c
...done. Time Elapsed: 104.218 ms.

   Created particles in .particles: 
   ParticleGroup with 100000 particles with total charge 1.0000000000000003e-11 C

In [6]:

P.plot('p')

P.plot('p')

In [7]:

P.plot('x', 'px')

P.plot('x', 'px')

In [8]:

P.plot('pz')

P.plot('pz')

In [9]:

P.plot('px')

P.plot('px')

In [10]:

P.plot('py')

P.plot('py')

In [11]:

from matplotlib import pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection='3d')

px, py, pz = P['px'], P['py'], P['pz']

p = np.sqrt(px**2 + py**2 + pz**2)

theta = np.arctan2(py, px)
phi = np.arccos(pz/p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

x = px/p
y = py/p
z = pz/p

ax.scatter(x[::200], y[::200], z[::200], '.');
ax.set_xlabel(r'$\hat{p}_x$')
ax.set_ylabel(r'$\hat{p}_y$')
ax.set_zlabel(r'$\hat{p}_z$')

from matplotlib import pyplot as plt fig = plt.figure() ax = fig.add_subplot(projection='3d') px, py, pz = P['px'], P['py'], P['pz'] p = np.sqrt(px**2 + py**2 + pz**2) theta = np.arctan2(py, px) phi = np.arccos(pz/p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 x = px/p y = py/p z = pz/p ax.scatter(x[::200], y[::200], z[::200], '.'); ax.set_xlabel(r'$\hat{p}_x$') ax.set_ylabel(r'$\hat{p}_y$') ax.set_zlabel(r'$\hat{p}_z$')

Out[11]:

Text(0.5, 0, '$\\hat{p}_z$')

In [12]:

theta = np.arctan2(py, px)
phi = np.arccos(pz/p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

plt.plot(pcs, hist);

theta = np.arctan2(py, px) phi = np.arccos(pz/p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 plt.plot(pcs, hist);

In [13]:

rns = np.random.random(100000)

rns = np.random.random(100000)

In [14]:

pa = 0
pb = np.pi

Ca = np.cos(pa)
Cb = np.cos(pb)

pa = 0 pb = np.pi Ca = np.cos(pa) Cb = np.cos(pb)

In [15]:

rho = 1/(Ca-Cb)

rho = 1/(Ca-Cb)

In [16]:

phis = np.linspace(pa, pb, 1000)

phis = np.linspace(pa, pb, 1000)

In [17]:

plt.plot(phis, rho*np.sin(phis));

plt.plot(phis, rho*np.sin(phis));

In [18]:

np.trapz(rho*np.sin(phis), phis)

np.trapz(rho*np.sin(phis), phis)

/tmp/ipykernel_2251/4275813628.py:1: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
  np.trapz(rho*np.sin(phis), phis)

Out[18]:

np.float64(0.999999175885426)

In [19]:

cdf = (Ca - np.cos(phis))*rho

cdf = (Ca - np.cos(phis))*rho

In [20]:

from scipy.integrate import cumulative_trapezoid as cumtrapz

from scipy.integrate import cumulative_trapezoid as cumtrapz

In [21]:

plt.plot(phis, cdf, phis, cumtrapz(rho*np.sin(phis), phis, initial=0));

plt.plot(phis, cdf, phis, cumtrapz(rho*np.sin(phis), phis, initial=0));

In [22]:

ps = np.arccos( Ca - rns*(Ca-Cb) )

ps = np.arccos( Ca - rns*(Ca-Cb) )

In [23]:

hist, pedges = np.histogram(ps, bins=30, density=True)

hist, pedges = np.histogram(ps, bins=30, density=True)

In [24]:

pcs = (pedges[1:] + pedges[:-1]) / 2

pcs = (pedges[1:] + pedges[:-1]) / 2

In [25]:

plt.plot(pcs, hist, phis, rho*np.sin(phis));

plt.plot(pcs, hist, phis, rho*np.sin(phis));

In [26]:

gen = Generator('data/beer.can.in.yaml', verbose=1)

#print(gen)

gen = Generator('data/beer.can.in.yaml', verbose=1) #print(gen)

In [27]:

P1 = gen.run()

P1 = gen.run()

Distribution format: gpt
Output file: beer.can.out.txt

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 10 pC.
   Number of macroparticles: 200000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 2 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   t distribution: uniform
      min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t:  1.1547 ps
   px distribution: Gaussian
      avg_px = 0 eV/c, sigma_px = 276.857 eV/c
   py distribution: Gaussian
      avg_py = 0 eV/c, sigma_py = 276.857 eV/c
   pz distribution: Gaussian
      avg_pz = 0 eV/c, sigma_pz = 276.857 eV/c
   Shifting avg_x = -2.10539E-05 mm -> 0 mm
   Scaling sigma_x = 0.999982 mm -> 1 mm
   Shifting avg_y = -1.97648E-06 mm -> 0 mm
   Scaling sigma_y = 0.999998 mm -> 1 mm
   Shifting avg_px = -0.0212516 eV/c -> 0 eV/c
   Scaling sigma_px = 276.849 eV/c -> 276.857 eV/c
   Shifting avg_py = -0.0259627 eV/c -> 0 eV/c
   Scaling sigma_py = 276.844 eV/c -> 276.857 eV/c
   Shifting avg_pz = -0.0372953 eV/c -> 0 eV/c
   Scaling sigma_pz = 276.843 eV/c -> 276.857 eV/c
   Shifting avg_t = -4.90022E-05 ps -> 0 ps
   Scaling sigma_t = 1.1547 ps -> 1.1547 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 220.904 eV/c, sigma_pz -> 166.887 eV/c
...done. Time Elapsed: 191.009 ms.

   Created particles in .particles: 
   ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C

In [28]:

P1.plot('x', 'px')

P1.plot('x', 'px')

In [29]:

P1.plot('p')

P1.plot('p')

In [30]:

P1.plot('pz')

P1.plot('pz')

In [31]:

P1['sigma_px'], P1['sigma_py'], P1['sigma_pz']

P1['sigma_px'], P1['sigma_py'], P1['sigma_pz']

Out[31]:

(np.float64(276.8570795554991),
 np.float64(276.8570795554991),
 np.float64(166.88683722575215))

In [32]:

gen = Generator('data/maxwell_boltzmann.beer.can.in.yaml', verbose=1)

gen = Generator('data/maxwell_boltzmann.beer.can.in.yaml', verbose=1)

In [33]:

P2 = gen.run()

P2 = gen.run()

Distribution format: gpt
Output file: beer.can.out.txt

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 10 pC.
   Number of macroparticles: 200000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 2 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   p distribution: Maxwell-Boltzmann
      p scale = 276.857 eV/c
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   t distribution: uniform
      min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t:  1.1547 ps
   Shifting avg_x = 1.63849E-05 mm -> 0 mm
   Scaling sigma_x = 0.999984 mm -> 1 mm
   Shifting avg_y = -2.77139E-06 mm -> 0 mm
   Scaling sigma_y = 0.999996 mm -> 1 mm
   Shifting avg_px = -0.000954618 eV/c -> -2.70524E-14 eV/c
   Scaling sigma_px = 276.853 eV/c -> 276.857 eV/rad⁰⋅⁵/c
   Shifting avg_py = 0.000856608 eV/c -> 0 eV/c
   Scaling sigma_py = 276.848 eV/c -> 276.857 eV/rad⁰⋅⁵/c
   Shifting avg_pz = -0.0142766 eV/c -> 2.70524E-14 eV/c
   Scaling sigma_pz = 276.833 eV/c -> 276.857 eV/rad⁰⋅⁵/c
   Shifting avg_t = -4.90022E-05 ps -> 0 ps
   Scaling sigma_t = 1.1547 ps -> 1.1547 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 220.909 eV/c, sigma_pz -> 166.881 eV/c
...done. Time Elapsed: 187.787 ms.

   Created particles in .particles: 
   ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C

In [34]:

P2['sigma_px'], P2['sigma_py'], P2['sigma_pz']

P2['sigma_px'], P2['sigma_py'], P2['sigma_pz']

Out[34]:

(np.float64(276.85707955549907),
 np.float64(276.85707955549907),
 np.float64(166.88083659251544))

In [35]:

gen = Generator('data/maxwell_boltzmann_KE.beer.can.in.yaml', verbose=1)

gen = Generator('data/maxwell_boltzmann_KE.beer.can.in.yaml', verbose=1)

In [36]:

P3 = gen.run()

P3 = gen.run()

Distribution format: gpt
Output file: beer.can.out.txt

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 10 pC.
   Number of macroparticles: 200000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 2 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   KE distribution: Maxwell-Boltzmann Energy
      kT = 150 meV

   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   t distribution: uniform
      min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t:  1.1547 ps
   Shifting avg_x = 1.63849E-05 mm -> 0 mm
   Scaling sigma_x = 0.999984 mm -> 1 mm
   Shifting avg_y = -2.77139E-06 mm -> 0 mm
   Scaling sigma_y = 0.999996 mm -> 1 mm
   Shifting avg_px = -1.75487E-09 meV·s/m -> -9.02058E-20 meV·s/m
   Scaling sigma_px = 0.000922998 meV·s/m -> 0.000922972 meV·s/m
   Shifting avg_py = 1.21968E-09 meV·s/m -> 0 meV·s/m
   Scaling sigma_py = 0.000922978 meV·s/m -> 0.000922972 meV·s/m
   Shifting avg_pz = -4.48933E-08 meV·s/m -> 9.02058E-20 meV·s/m
   Scaling sigma_pz = 0.00092294 meV·s/m -> 0.000922972 meV·s/m
   Shifting avg_t = -4.90022E-05 ps -> 0 ps
   Scaling sigma_t = 1.1547 ps -> 1.1547 ps
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 0.0007366 meV·s/m, sigma_pz -> 0.000556145 meV·s/m
...done. Time Elapsed: 190.986 ms.

   Created particles in .particles: 
   ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C

In [37]:

P3.plot('kinetic_energy')

P3.plot('kinetic_energy')

In [38]:

P3.plot('x', 'px')

P3.plot('x', 'px')

In [39]:

P3['sigma_px'], P3['sigma_py'], P3['sigma_pz']

P3['sigma_px'], P3['sigma_py'], P3['sigma_pz']

Out[39]:

(np.float64(276.70000872169317),
 np.float64(276.70000872169317),
 np.float64(166.7279863589753))

In [40]:

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

px, py, pz = P3['px'], P3['py'], P3['pz']

p = np.sqrt(px**2 + py**2 + pz**2)

theta = np.arctan2(py, px)
phi = np.arccos(pz/p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

x = px/p
y = py/p
z = pz/p

ax.scatter(x[::500], y[::500], z[::500], '.');
ax.set_xlabel(r'$\hat{p}_x$')
ax.set_ylabel(r'$\hat{p}_y$')
ax.set_zlabel(r'$\hat{p}_z$')

fig = plt.figure() ax = fig.add_subplot(projection='3d') px, py, pz = P3['px'], P3['py'], P3['pz'] p = np.sqrt(px**2 + py**2 + pz**2) theta = np.arctan2(py, px) phi = np.arccos(pz/p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 x = px/p y = py/p z = pz/p ax.scatter(x[::500], y[::500], z[::500], '.'); ax.set_xlabel(r'$\hat{p}_x$') ax.set_ylabel(r'$\hat{p}_y$') ax.set_zlabel(r'$\hat{p}_z$')

Out[40]:

Text(0.5, 0, '$\\hat{p}_z$')

In [41]:

KE = np.linspace(0, 150, 10000)

E0=0

A1 = 0.8,
m1 = 8

A2 = 0.1
m2 = 90

PKE = A1*np.exp(-np.abs(KE-E0)/m1) + A2*np.exp(-np.abs(KE-E0)/m2) 
PKE = PKE/np.trapz(PKE, KE)  # Numerically intergate to normalize

plt.plot(KE, PKE);
plt.xlabel('KE (eV)');
plt.ylabel('$\\rho(KE)$ (1/eV)');

KE = np.linspace(0, 150, 10000) E0=0 A1 = 0.8, m1 = 8 A2 = 0.1 m2 = 90 PKE = A1*np.exp(-np.abs(KE-E0)/m1) + A2*np.exp(-np.abs(KE-E0)/m2) PKE = PKE/np.trapz(PKE, KE) # Numerically intergate to normalize plt.plot(KE, PKE); plt.xlabel('KE (eV)'); plt.ylabel('$\\rho(KE)$ (1/eV)');

/tmp/ipykernel_2251/1250638860.py:12: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
  PKE = PKE/np.trapz(PKE, KE)  # Numerically intergate to normalize

In [42]:

dat = np.zeros((len(KE),2))
dat[:,0], dat[:,1] = KE, PKE

np.savetxt('KEdist.txt', dat, header="KE     PKE", comments='')

dat = np.zeros((len(KE),2)) dat[:,0], dat[:,1] = KE, PKE np.savetxt('KEdist.txt', dat, header="KE PKE", comments='')

In [43]:

q = 1000*1.60217663e-19

print(q)

input_yaml = """
n_particle: 1000
species: electron
start: 
  type: cathode
random:
  type: hammersley
total_charge:
  units: C
  value: 1.60217663e-17
r_dist:
  max_r:
    units: mm
    value: 14.6
  type: radial_uniform
KE_dist:
  file: KEdist.txt
  units: eV
  type: file1d
transforms:
  sx:
    avg_x:
      units: millimeter
      value: -50
    type: set_avg x
  sy:
    avg_y:
      units: millimeter
      value: 50
    type: set_avg y
  sz:
    avg_z:
      units: millimeter
      value: 65
    type: set_avg z
output:
    file: None
"""

q = 1000*1.60217663e-19 print(q) input_yaml = """ n_particle: 1000 species: electron start: type: cathode random: type: hammersley total_charge: units: C value: 1.60217663e-17 r_dist: max_r: units: mm value: 14.6 type: radial_uniform KE_dist: file: KEdist.txt units: eV type: file1d transforms: sx: avg_x: units: millimeter value: -50 type: set_avg x sy: avg_y: units: millimeter value: 50 type: set_avg y sz: avg_z: units: millimeter value: 65 type: set_avg z output: file: None """

1.6021766299999998e-16

In [44]:

import yaml

inputs = yaml.safe_load(input_yaml)

n = 1000

inputs['n_particle'] = n
inputs['output']['file'] = 'test_ion_writer.ion'
inputs['output']['type'] = 'simion'
inputs['total_charge']['value'] = n*1.60217663e-19
inputs['transforms']['sx']['avg_x']['value'] = -50
inputs['transforms']['sy']['avg_y']['value'] = 50
inputs['transforms']['sz']['avg_z']['value'] = 55

import yaml inputs = yaml.safe_load(input_yaml) n = 1000 inputs['n_particle'] = n inputs['output']['file'] = 'test_ion_writer.ion' inputs['output']['type'] = 'simion' inputs['total_charge']['value'] = n*1.60217663e-19 inputs['transforms']['sx']['avg_x']['value'] = -50 inputs['transforms']['sy']['avg_y']['value'] = 50 inputs['transforms']['sz']['avg_z']['value'] = 55

In [45]:

gen = Generator(inputs, verbose=1)

gen = Generator(inputs, verbose=1)

In [46]:

P = gen.run()

P = gen.run()

Distribution format: simion
Output file: /home/runner/work/distgen/distgen/docs/examples/test_ion_writer.ion

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 1.60218E-16 C.
   Number of macroparticles: 1000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 14.6 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   KE distribution: KE-distribution file: "/home/runner/work/distgen/distgen/docs/examples/KEdist.txt"
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   Shifting avg_x = 0.00564031 mm -> 0 mm
   Scaling sigma_x = 7.29551 mm -> 7.3 mm
   Shifting avg_y = 0.00120927 mm -> 0 mm
   Scaling sigma_y = 7.29741 mm -> 7.3 mm
   Shifting avg_px = -6.09774E-09 eV·s/m -> -9.89015E-22 eV·s/m
   Scaling sigma_px = 1.11822E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_py = -2.50205E-09 eV·s/m -> 0 eV·s/m
   Scaling sigma_py = 1.11987E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_pz = -2.52855E-08 eV·s/m -> 9.89015E-22 eV·s/m
   Scaling sigma_pz = 1.11509E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 8.07962E-06 eV·s/m, sigma_pz -> 7.72348E-06 eV·s/m
   Applying user supplied transform: "sx" = set_avg x...
      Setting avg_x -> -50 mm.
   Applying user supplied transform: "sy" = set_avg y...
      Setting avg_y -> 50 mm.
   Applying user supplied transform: "sz" = set_avg z...
      Setting avg_z -> 55 mm.
...done. Time Elapsed: 26.2625 ms.

   Created particles in .particles: 
   ParticleGroup with 1000 particles with total charge 1.6021766300000003e-16 C

In [47]:

P.plot('x', 'y', figsize=(5,5))

P.plot('x', 'y', figsize=(5,5))

In [48]:

P.plot('x', 'px', figsize=(5,5))

P.plot('x', 'px', figsize=(5,5))

In [49]:

P.plot('kinetic_energy')

P.plot('kinetic_energy')

In [50]:

from distgen.writers import writer

writer('simion', gen.beam(), 'test_ion_writer.ion', params={'color':0})

from distgen.writers import writer writer('simion', gen.beam(), 'test_ion_writer.ion', params={'color':0})

Distribution format: simion
Output file: /home/runner/work/distgen/distgen/docs/examples/test_ion_writer.ion

Creating beam distribution....
   Beam starting from: cathode
   Total charge: 1.60218E-16 C.
   Number of macroparticles: 1000.
   Assuming cylindrical symmetry...
   r distribution: radial uniform
      min_r = 0 mm, max_r = 14.6 mm
   theta distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   KE distribution: KE-distribution file: "/home/runner/work/distgen/distgen/docs/examples/KEdist.txt"
   azimuthal angle distribution: uniform theta
      min_theta = 0 rad, max_theta = 6.28319 rad
   polar angle distribution: uniform phi
      min_phi = 0 rad, max_phi = 3.14159 rad
   Shifting avg_x = 0.00564031 mm -> 0 mm
   Scaling sigma_x = 7.29551 mm -> 7.3 mm
   Shifting avg_y = 0.00120927 mm -> 0 mm
   Scaling sigma_y = 7.29741 mm -> 7.3 mm
   Shifting avg_px = -6.09774E-09 eV·s/m -> -9.89015E-22 eV·s/m
   Scaling sigma_px = 1.11822E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_py = -2.50205E-09 eV·s/m -> 0 eV·s/m
   Scaling sigma_py = 1.11987E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Shifting avg_pz = -2.52855E-08 eV·s/m -> 9.89015E-22 eV·s/m
   Scaling sigma_pz = 1.11509E-05 eV·s/m -> 1.11773E-05 eV·s/m
   Cathode start: fixing pz momenta to forward hemisphere
      avg_pz -> 8.07962E-06 eV·s/m, sigma_pz -> 7.72348E-06 eV·s/m
   Applying user supplied transform: "sx" = set_avg x...
      Setting avg_x -> -50 mm.
   Applying user supplied transform: "sy" = set_avg y...
      Setting avg_y -> 50 mm.
   Applying user supplied transform: "sz" = set_avg z...
      Setting avg_z -> 55 mm.
...done. Time Elapsed: 26.4134 ms.

In [51]:

from scipy.constants import physical_constants

mc2 = 1e6 * physical_constants['electron mass energy equivalent in MeV'][0]
e_= physical_constants['elementary charge'][0]
me = physical_constants['electron mass in u'][0]

def particle_group_to_SIMION(P, filename, color=0):
    
    header=';0'
    
    simion_params= ['TOB', 'MASS', 'CHARGE', 'X', 'Y', 'Z', 'AZ', 'EL', 'KE', 'CWF', 'COLOR']
    
    simion_units = {'TOB':'usec', 'MASS':'amu', 'CHARGE':'e', 'X':'mm', 'Y':'mm', 'Z':'mm', 'AZ':'deg', 'EL':'deg', 'CWF':'', 'COLOR':''}
    
    data = np.zeros( (len(P), len(simion_params)) )
    
    data[:, simion_params.index('TOB')] = P.t*1e6    # [P.t] = sec, convert to usec
    
    if(P.species == 'electron'):
        data[:, simion_params.index('MASS')] = np.full(len(P), me)
        data[:, simion_params.index('CHARGE')] = np.full(len(P), -1)
    else:
        raise ValueError(f'Species {P.species} is not supported')
    
    data[:, simion_params.index('X')] =  P.z*1e3
    data[:, simion_params.index('Y')] =  P.y*1e3
    data[:, simion_params.index('Z')] = -P.x*1e3
    
    px =  P.pz
    py =  P.py
    pz = -P.px
    
    data[:, simion_params.index('KE')] = P.kinetic_energy       # [eV] 
    data[:, simion_params.index('AZ')] = np.arctan2(-pz, px) * (180/np.pi)   # [deg]
    data[:, simion_params.index('EL')] = np.arctan2(py, np.sqrt(px**2 + pz**2) ) * (180/np.pi) # [deg]
    
    data[:, simion_params.index('CWF')] = P.weight/e_                       # Charge Weighting Factor, derive from particle group weights
    data[:, simion_params.index('COLOR')] = np.full(len(P), color)
    #fname, X, fmt='%.18e', delimiter=' '
    
    np.savetxt(filename, data, delimiter=',', header=header, comments='', fmt='  %.9e')

from scipy.constants import physical_constants mc2 = 1e6 * physical_constants['electron mass energy equivalent in MeV'][0] e_= physical_constants['elementary charge'][0] me = physical_constants['electron mass in u'][0] def particle_group_to_SIMION(P, filename, color=0): header=';0' simion_params= ['TOB', 'MASS', 'CHARGE', 'X', 'Y', 'Z', 'AZ', 'EL', 'KE', 'CWF', 'COLOR'] simion_units = {'TOB':'usec', 'MASS':'amu', 'CHARGE':'e', 'X':'mm', 'Y':'mm', 'Z':'mm', 'AZ':'deg', 'EL':'deg', 'CWF':'', 'COLOR':''} data = np.zeros( (len(P), len(simion_params)) ) data[:, simion_params.index('TOB')] = P.t*1e6 # [P.t] = sec, convert to usec if(P.species == 'electron'): data[:, simion_params.index('MASS')] = np.full(len(P), me) data[:, simion_params.index('CHARGE')] = np.full(len(P), -1) else: raise ValueError(f'Species {P.species} is not supported') data[:, simion_params.index('X')] = P.z*1e3 data[:, simion_params.index('Y')] = P.y*1e3 data[:, simion_params.index('Z')] = -P.x*1e3 px = P.pz py = P.py pz = -P.px data[:, simion_params.index('KE')] = P.kinetic_energy # [eV] data[:, simion_params.index('AZ')] = np.arctan2(-pz, px) * (180/np.pi) # [deg] data[:, simion_params.index('EL')] = np.arctan2(py, np.sqrt(px**2 + pz**2) ) * (180/np.pi) # [deg] data[:, simion_params.index('CWF')] = P.weight/e_ # Charge Weighting Factor, derive from particle group weights data[:, simion_params.index('COLOR')] = np.full(len(P), color) #fname, X, fmt='%.18e', delimiter=' ' np.savetxt(filename, data, delimiter=',', header=header, comments='', fmt=' %.9e')

In [52]:

particle_group_to_SIMION(P, 'text_ion_file.ion')

particle_group_to_SIMION(P, 'text_ion_file.ion')

In [53]:

def read_simion_ION_file(filename):
    
    data = np.loadtxt(filename, comments=';', delimiter=',', skiprows=1)
    
    simion_params= ['TOB', 'MASS', 'CHARGE', 'X', 'Y', 'Z', 'AZ', 'EL', 'KE', 'CWF', 'COLOR']
    
    return {simion_params[ii]:data[:,ii] for ii, p in enumerate(simion_params)}

def read_simion_ION_file(filename): data = np.loadtxt(filename, comments=';', delimiter=',', skiprows=1) simion_params= ['TOB', 'MASS', 'CHARGE', 'X', 'Y', 'Z', 'AZ', 'EL', 'KE', 'CWF', 'COLOR'] return {simion_params[ii]:data[:,ii] for ii, p in enumerate(simion_params)}

In [54]:

ions1 = read_simion_ION_file('text_ion_file.ion')
ions2 = read_simion_ION_file('test_ion_writer.ion')

ions1 = read_simion_ION_file('text_ion_file.ion') ions2 = read_simion_ION_file('test_ion_writer.ion')

In [55]:

fig = plt.figure()
ax = fig.add_subplot(projection='3d')

px, py, pz = P['px'], P['py'], P['pz']

p = np.sqrt(px**2 + py**2 + pz**2)

theta = np.arctan2(py, px)
phi = np.arccos(pz/p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2

x = px/p
y = py/p
z = pz/p

ax.scatter(x[::], y[::], z[::], '.');
ax.set_xlabel(r'$\hat{p}_x$')
ax.set_ylabel(r'$\hat{p}_y$')
ax.set_zlabel(r'$\hat{p}_z$')

fig = plt.figure() ax = fig.add_subplot(projection='3d') px, py, pz = P['px'], P['py'], P['pz'] p = np.sqrt(px**2 + py**2 + pz**2) theta = np.arctan2(py, px) phi = np.arccos(pz/p) hist, pedges = np.histogram(phi, bins=100, density=True) pcs = (pedges[1:] + pedges[:-1]) / 2 x = px/p y = py/p z = pz/p ax.scatter(x[::], y[::], z[::], '.'); ax.set_xlabel(r'$\hat{p}_x$') ax.set_ylabel(r'$\hat{p}_y$') ax.set_zlabel(r'$\hat{p}_z$')

Out[55]:

Text(0.5, 0, '$\\hat{p}_z$')

In [56]:

for p in ions1: print(f'{p}:', max(np.abs(ions1[p] - ions2[p])))

for p in ions1: print(f'{p}:', max(np.abs(ions1[p] - ions2[p])))

TOB: 0.0
MASS: 0.0
CHARGE: 0.0
X: 0.0
Y: 0.0
Z: 0.0
AZ: 0.0
EL: 0.0
KE: 0.0
CWF: 0.0
COLOR: 0.0

In [57]:

plt.plot(ions1['KE'], ions1['KE']/ions2['KE'], '.');
plt.xlabel('KE (eV)');

plt.plot(ions1['KE'], ions1['KE']/ions2['KE'], '.'); plt.xlabel('KE (eV)');

In [58]:

os.remove('text_ion_file.ion')
os.remove('KEdist.txt')
os.remove('test_ion_writer.ion')

os.remove('text_ion_file.ion') os.remove('KEdist.txt') os.remove('test_ion_writer.ion')

In [ ]:

In [ ]: